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Broken Symmetry in Ideal Magnetohydrodynamic Turbulence 


John V. Shebalin* 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23681, USA 

A numerical study of the long-time evolution of a number of cases of inviscid, isotropic, incompressible, 
three-dimensional fluid and magneto-fluid turbulence has been completed. The results confirm that ideal 
magnetohydrodynamic turbulence is non-ergodic if there is no external magnetic field present. This is due 
essentially to a canonical symmetry being broken in an arbitrary dynamical representation. The broken 
symmetry manifests itself as a coherent structure, t. e., a non-zero time-averaged part of the turbulent 
magnetic field. The coherent structure is observed, in one case, to contain about eighteen percent of the 
total energy. 


* Research supported by the National Aeronautics and Space Administration. This work was performed while the author 
was in residence as a Visiting Scientist at the Institute for Computer Applications in Science and Engineering (ICASE), NASA 
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1 Introduction 


Several years ago, a numerical study of isotropic, inviscid, incompressible, two-dimensional (2D) turbu- 
lence indicated that ideal magnetohydrodynamic (MHD) turbulence was non-ergodic [1]. A Fourier spectral 
method was used, and the non-ergodicity manifested itself in the appearance of significant, non-zero time- 
averages for the Fourier coefficients of the turbulent magnetic field. These results contrasted with the 
predictions of canonical ensemble theory [2] and indicated the presence of a dynamically broken symmetry. 

Here, this work is extended in two ways. First, the dynamics become fully three-dimensional (3D), 
and second, a more accurate (third-order) time-integration scheme is used (a comparison is made with the 
second-order scheme used previously). The results are qualitatively the same as previously seen, although a 
stronger effect is observed. Also, the source of the broken symmetry is more explicitly identified. 

2 Basic Equations 

The equations which describe ideal, incompressible 3D fluid and magneto-fluid dynamics are 

d t B = V x (u x B) (1) 

= Vx(uxw+jxB) (2) 

Here, the fluid velocity is u, the vorticity is w = V x u, the magnetic field is B = b + B 0 (where b is 
the turbulent part and B 0 is a constant, externally imposed part), and the current is j = V x b. Also, 
V • u = V • B = 0. 

If B = 0 for all time, then (I) is identically satisfied and (2) assumes the vorticity form of the incom- 
pressible Euler equation. If B<? = 0 while b ^ 0, then (1) and (2) are unchanged, i.e., symmetric, under the 

substitution b — ► — b. If, however, B 0 ^ 0 also, then the equations are no longer symmetric under b — ► — b. 

Note that these equations are inviscid. It is well known that solutions to these equations are quite 
different from those for dissipative equations, even if the viscosity and resistivity become very small, so long 
as they are not identically zero. However, it is possible to determine a priori statistical solutions in the 
inviscid case (as will be discussed presently), while it has not proven possible to do this in the dissipative 
case. 
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This provides for a very practical use of the results presented herein, that is, as a test of the correctness 
of any computer code which attempts to simulate incompressible, homogeneous, dissipative turbulence. In 
any such code, the viscosity (and resistivity) can be set identically to zero; the statistical properties of a 
numerical solution must then be described by the inviscid theory described here. There are, however, some 
refinements which must be incorporated into existing theory; these refinements, the existence of broken 
symmetry and non-ergodicity in ideal (3D) MHD turbulence, are the subject of this paper. 


3 Canonical Ensembles 


The theory of canonical or ‘absolute equilibrium 5 ensembles, as it applies to homogeneous turbulence, has 
been discussed extensively before [1, 3] (and the many references therein). In brief, the physical variables u, 
b, a?, and j are expanded in truncated Fourier series, for example: 


■i & mo x 

b ( x >*) = b (MK k : 


( 3 ) 


Here, the sum (and each similar sum appearing henceforth) is over all k such that k m in < |k| < k ma x < N/2, 
where N is the number of points in each of the three spatial dimensions. Since the various fields, such as b(x), 
are real, their coefficients satisfy b(k) = b*(-k), where £ *’ denotes complex conjugation (here and henceforth, 
explicit time-dependence will be dropped for brevity). Also, u(k) = ik ~ 2 k x u?(k) and j(k) = ik x b(k), 
and V * b(x) = 0 — ► ik * b(k) = 0, etc. 

The independent real and imaginary parts of the coefficients u(k) and b(k) can be used to label the axes 
of a multidimensional phase space. The corresponding dynamical system is described by a single point in 
the phase space, which moves about as the system evolves in time. The probability that the system point 
is in any part of phase space can be described by a canonical distribution function which depends only on a 
small set of conserved quantities, the integral invariants of the dynamical system. Once the joint probability 
distribution has been found, then the equilibrium energy spectra (kinetic and magnetic) are determined, 
even though u(k) and b(k) are random variables. 

In the case of isotropic, incompressible 3D Euler turbulence (B = 0), the integral invariants are the 
energy E and the kinetic helicity H k . For ideal, isotropic, incompressible 3D MHD turbulence with B c = 0, 
the integral invariants are the energy E } the cross helicity H c , and the magnetic helicity H m \ if B 0 / 0, then 
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H m is no longer conserved, although E and H c are still integral invariants. These three different situations 
will now be presented in more detail. 


3.1 Euler Turbulence 


The integral invariants for incompressible 3D Euler turbulence [4, 5] are 

1 t m a x -i ^ max 

B = jjk E MOI’ ” d H ‘=2N* E“< k >'"'( k >' 

fcmin 

while the distribution function is 

/ a 2 _ d' 2 k 2 \ ^ 

D = Cexp{-aE-0H k ) where C = JJ ( • 


(4) 


(5) 


The product above (and each product appearing henceforth) is taken over all k such that k m i n < |k| < k max > 
but only over independent values of k, i.e., if k is used, then — k is not. Using this distribution function, an 
expectation value can be defined as 


max 

QD n d 3 u(k), 

JU 

*■ m i n 


( 6 ) 


with which expectation values of the moments of the Fourier coefficients can be found: 
(u*(k)> = (u/(k)) = 0, (|ufl(k)| 2 > = <| U/ (k)| 2 > = 2(aa _-^ -py 

— ^N 3 Bk 2 

and <u fi (k) • w fl (k)} = (u/(k) • w/(k)> = 2 ^ q2 _^ 


(7) 


where the subscripts R and I denote real and imaginary parts, respectively, of the complex coefficients, and 
where all the components of u(k) for a given k have equal expectation values. A straightforward algebraic 
manipulation of these expectation values yields 

Z(Q) 


a — 


(E) (Q) - (H k ) 


a a ( H *> 

2 and /? = _ w a ’ 


( 8 ) 


where 


* fc max i ^ mar 

z =Em1L l and fc2 l u ( k )l 2 - 


27V 3 ^ ‘ ” 27V 3 

min ^min 

Note that the expected value of the mean squared vorticity (the enstrophy ft) is determined once the 
expectation values (7) of the |u(k)| 2 are known. 
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Although the enstrophy Q is not a conserved quantity, it does have an expectation value, so that the 
parameters a and 0 are functions of only one unknown quantity, (fi). The quantity Z is the ratio of the 
number of independent values of k to iV 3 , the total number of points in the spatial grid. 


3.2 Ideal MHD Turbulence 

The integral invariants of ideal, incompressible 3D MHD turbulence with B 0 = 0 [6] are 

■i kmax -i ^mn 

E = 2^3 £ [Kk)! ! + IMk)| ! l . He = 2^3 £ »<*) ■ b '(k). 

k m in ^min 


and 


-i * max 


^m»n 


while the distribution function is 


D = Cexp (— otE — 0H C — 7i7 m ), where C = 

^mtn L 


(a 2 - /? 2 /4) 2 - a 2 j 2 fk 2 




t r 4 iV 12 


(9) 


( 10 ) 


(Please note that the formulas appearing in this subsection correct those which appeared previously [1].) 

Using this distribution function, expectation values (where (6) now includes integration over the inde- 
pendent components of the b(k)) for moments of the Fourier coefficients can be found: 


(Ufl(k)) = <u/(k)> = 0, <|u fi (k)| 2 ) = (|u,(k)| 2 ) = - N 3 


a (a 2 — /? 2 /4 — ~f 2 /k 2 ) 


[(a 2 — /? 2 /4) —a 2 j 2 /k 2 j 


(bfl(k)} = (b,(k)> = 0, (|b R (k)| 2 ) = <|b/(k)| 2 ) = -N 3 


a (a 2 — /? 2 /4) 


(a 2 — P 2 /4 y — a 2 j 2 /k 2 


{u fl (k)-bj l (k)) = (u / (k)-bj(k ))=^N 3 


-p (a 2 -/? 2 / 4) 


(a 2 — P 2 /4y — a 2 ~t 2 /k 2 


an(k).bi(k)> = (j/(k)-bj(k)> = ^ 3 


o 

—a 7 


(ii) 


2 [(a 2 -p 2 /4f-a 2 'r 2 /k 2 \’ 

where all the components of u(k) for a given k have equal expectation values, and where all the components 
of b(k) for a given k also have equal expectation values. Note that the energy spectra in the MHD case (11) 
peak at low k = |k|, while in the Euler case (7), it peaks at high k . 
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A straightforward algebraic manipulation of these expectation values yields 


a = 


ZR(R+1){E) 


2 R(E) 2 -(R + iy {H c y 


0 =- 2 


{R+l){H e ) 
R(E) ' 
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where 


{Em 

(£*}' 


( 12 ) 


Thus a, /?, and 7 are functions of only one unknown quantity, iZ, which is the ratio of the expectation values 
of the magnetic part of the total energy to kinetic part of the total energy. 

The integral invariants of ideal, incompressible 3D MHD turbulence with B 0 ^ 0 are only E and H Ci as 
given above, while H m is no longer conserved. In this case, R = 1 and 7 = 0. 


3.3 Distribution Parameters 

The parameters a, /?, and 7, which appear in the canonical distribution functions, must be given specific 
values in order that definite predictions can be made. This can be done in a number of ways. First, it 
can be assumed that the integral invariants of the various cases have fixed values (although they actually 
fluctuate slightly in the time evolution of a canonical system). The distribution parameters then vary only 
with respect to a single quantity: (ft) for 3D Euler turbulence and R for ideal 3D MHD turbulence (except 
for the case where B 0 ^ 0, where R = 1). The equilibrium entropy [1], 5 = S 0 — logC, where C is the 
normalizing coefficient of the distribution function, is thus also a function of a single quantity. Since S’ is a 
minimum with respect to the distribution parameters [7], then this fact can be used to uniquely determine 
(ft) or R, as the case may be. Note that this can be done prior to any numerical simulation of the time 
evolution of an isotropic, inviscid turbulent system. 

A second method is to run a numerical simulation and actually determine the time-averaged values Q 
of all the quantities Q which are required, so that Q can be substituted for (Q). In particular, ft — ► (ft) or 
R R. A third, and more practical method, is to use the time-averaged values of the necessary integral 
invariants, but determine (ft) or R through a least-squares procedure. This third method will turn out to 
be the most useful. 

Once the canonical predictions of the moments of u(k) and b(k) are available, and numerical simulations 
have produced corresponding time- aver ages, then comparisons can be made. This is essentially a test of the 
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ergodicity of the various dynamical systems, z.e., the equivalence of ensemble and time averages. Numerical 
results, including these time averages, are presented next, followed by a discussion addressing any anomalies 
which are observed. 

4 Numerical Results 

As in the previous work [1], a Fourier spectral method [9] with shifted-grid dealiasing [10] was used to 
solve the dynamical equations numerically. In the previous effort, however, the time-integration method 
was a second-order Runge-Kutta scheme (RK2), while the method primarily employed in this work was a 
third-order scheme, consisting of an Adams-Bashforth predictor coupled with an Adams-Moulton corrector 
(AB3) [8]. The advantages of AB3 over RK*2 are higher accuracy and higher speed, although slightly more 
storage was required in computer memory. In order to better relate the 3D results obtained here to the 2D 
ones obtained previously [1], two 3D runs which had used AB3 were partially re-run using RK2, with initial 
conditions remaining the same. 

The numerical simulations performed in the course of this work are presented in Table 1. All runs were 
done on a Cray YMP, using a 16 3 grid with = 56 and At = 0.001. Each run is designated by three 
characters: the first is either E (for Euler) or M (for MHD), the second is a number representing a different 
set of initial conditions, and the last denotes the time-integration method used, A (for AB3) or R (for RR2). 
(The ‘initial conditions 5 include the values of B 0 and k min , as well as {u>(k), b(k); k m i n < |k| < fc ma x} at 
t = 0.) Each of the integral invariants listed in Table 1 fluctuated no more than a few parts per million 
during their respective simulations. 

In each of the simulations, the initial spectra satisfied |u(k)| 2 ~ |b(k)| 2 ~ fc 4 exp(-2 k 2 fk 2 0 ) with k 0 = 2 
(although the exact shape of the initial spectra was not critical; the critical initial quantities were the values 
of the total energy and helicities). Time- averages of the components of u>(k) and b(k) as well as of their 
squares (|6 x /*(k)j 2 , R: real, etc.) were taken, averaging every 50 time steps. (There are six ‘components’ 
of w(k) and six of b(k), since the x— y y— , and ^-components of these each have a real and an imaginary 
part.) The second-order moments which correspond to the same value of k = |k| can be combined to produce 
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<1 


average modal energies 


^) = iEK k )i 2 and ^* )= db? |b(k)|2 ’ (i3) 

where the sum Yk on ^ over ^ same va I ue of k = |k| and where ra* is the number of wave 

vectors k in such a sum. The modal values of enstrophy and mean squared current were then fi = k 2 Ek(k) 
and J = k 2 E m (k) 1 respectively. 

These modal energies can be compared with those obtained from ensemble- averaging, once the various 
distribution parameters are found, as described previously. Remember that the distribution parameters for 
the Euler cases depended on the average value of the enstrophy f2, while those for the MHD cases depended 
on the average ratio of magnetic to kinetic energy R. For example, in the MHD run MIA, a time-average 
yields R t = 1.50795, equilibrium entropy minimization gives R e = 1.50799, and minimizing the root-mean- 
square (rms) error between time- and ensemble- aver age modal energy spectra produces i? rmj = 1.50599. 
Using the rms values in all cases gives the distribution parameters presented in Table 2. 

To see how well the time- and ensemble- averaged modal energy spectra compare, consider Figure 1, where 
representative spectra are shown (log = logio in the Figures). Figure la corresponds to case M2 A where 
k min = 1; the goodness of fit is typical of all the Euler and MHD runs. Figure lb corresponds to case M3 A 
where fc min = 2; the initial values of the u>(k) and b(k) for M3 A were the same as those of M2 A except that 
all of the |k| = 1 coefficients were set to zero and not allowed to grow during the simulation. These graphs 
show the average energy for modes with the same k = |k|; since there are twice as many modes with k = 2 
as with k = 1, Figures la and lb indicate that the modes with k = contain the same total energy. The 
exact value of k m i n is thus not qualitatively important in any given simulation. 

It is informative to consider the exact values taken by the Fourier modes during their time evolution. 
For example, the values of u>(k) and b(k) for k =(1,0,0) for run MIA are shown in Figure 2, where the time 
evolution is shown by plotting the real vs the imaginary parts of the coefficients for all times between t = 0 
and t = 1000. The behavior is decidedly non-ergodic since the components clearly do not have zero mean 
values, and thus do not match the ensemble prediction. It is interesting to note that although average modal 
behavior (Figure 1) is essentially as predicted, the detailed evolution of the modal coefficients is not. Also, 
notice in Figure 2 that we appear to have 

Wy ~ iw z and b y ~ ib z . (14) 
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This relationship was also apparent for all the other coefficients which were examined, as long as the coefficient 
mean values were sufficiently larger than the corresponding standard deviation. The relationship (14) is, in 
fact, b ~ V x b, i.e.Jxb^O. In other words, the ideal MHD has become force-free [11], apparently on a 
mode- by- mode basis. 

The average coefficients u?(k) and b(k) can be thought of as the coherent structure underlying these 
simulations of ideal, isotropic turbulence. These coefficients can be used to calculate the coherent energy 
present in the turbulence: 

2^ £l Q ( k )| 2 and £ mW=2^El 6 ( k )i 2 > ( 15 ) 

k k 

(The various helicities can also be calculated in a similar manner.) The average coefficients exist for all cases 
at t = 250, so that a comparison can be made; the results of this are presented in Table 3. Notice that the 
coherent quantities in Table 3 are a sizeable fraction of the corresponding total quantities in Table 1; for 
example, the coherent energy in run MIA is 18.8% of the total energy, while the coherent magnetic helicity 
in the same run is 83% of the total! 

For those runs which went beyond t = 250, the coherent energies between t = 520 and t = 750 are 
shown in Figure 3. It is evident that the MHD simulations with B 0 = 0 (MIA and M2A) have the largest 
coherent energies, and that this coherent energy is primarily magnetic. The MHD case with B 0 = 1 (M4A) 
has effectively no coherent energy, while in the Euler case (E2A), the coherent energy is clearly decaying 
(although still 3% of the total at f = 750). 

5 Broken Symmetry 

In a Fourier representation (3), the basic equations (1) and (2) take the form 

b(k) = ikx E Kp)xb(q)] +*(k-B„)u(k) (!6) 

p+q=k 

u>(k) = ikx E [u(p) x u»(q) + j(p) x b(q)j + i(k • B 0 ) j(k) (17) 

p+q=k 

These equations have the symmetry (i.e., invariance under) w(k) — ► e’ k a w(k), b(k) — ► e ,k a b(k) which 

merely reflects the isotropy of space: i.e., place x -+ x + a into (3). However, an average over the displace- 

ments a merely reproduces the known fact that the k = 0 modes of u and b are zero and remain so, i.e., do 
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not participate in the dynamics of the ideal MHD turbulence. 

In the case where B 0 = 0, the equations (16) again also have the symmetry b — ► -b; the import of this 
is that all members of the ensemble of possible realizations of ideal MHD turbulence are, in fact, paired: if 
(u(k), b(k); k m i n < |k| < k max } is a possible realization, then so is {u(k), —b(k); k m i n < W ^ kmax }• 

Thus what is calculated in the ensemble average (b(k)) is really ([b+(k) + b_(k)]/2), where b+(k) = 
b(k) and b_(k) = — b(k), in which case (b(k)) = 0 automatically. However, (b+(k)) and (b_(k)) are 
undefined; if they are non-zero, then the ensemble symmetry is said to be broken in a given realization (and 
the dynamics are non-ergodic). That they can be non-zero, and substantially so, has been demonstrated in 
the numerical results presented herein. 

6 Conclusion 

In this paper, a numerical study of 3D Euler and ideal MHD turbulence was presented. A number of 
different cases were presented in which the dynamic evolution was allowed to proceed for a considerable 
length of dimensionless simulation time. These results clearly demonstrate the presence of broken symmetry, 
non-ergodicity, and coherent structure in ideal 3D MHD turbulence. 

Realistic MHD turbulence (except in superfluids) is not ideal since resistive and viscous processes dissipate 
energy. The results presented here then address, rather than practical cases, only numerical simulations in 
which fluid and magnetic dissipation can be set identically equal to zero. For these simulations, however, 
the ideal theory provides a very useful means of testing computer codes. 

An interesting, open question is: Does the non-ergodic behavior seen in the ideal magnetohydrodynamic 
case also occur in real geophysical and astrophysical systems where dynamo activity has been observed [11]? 
An answer to this question, however, requires the introduction of dissipation into numerical simulations, 
along with a substantial increase in grid size. Although this was beyond the scope of the present work, it is 
a straightforward extension of it, provided that sufficiently large computers can be utilized. 

Additionally, it should be mentioned that there has been some very interesting work in dynamical systems 
theory concerning broken symmetry [12]. A full discussion of any connection to the numerical work presented 
here is also beyond the scope of the present work. 

Finally, I would like to thank Dr. M. Y. Hussaini for the opportunity to spend a year in residence at 
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ICASE as a Visiting Scientist, and to thank Professor Geoffrey Lilley for his useful comments concerning 
this paper. 
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Table 3. Coherency at t = 250 
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— 


— 

— 
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m 
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— 

K 
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— 



M3A 

HHH 
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0.0368 

M4A 


0.00148 

— 
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Figure 2. Components of two vector modes of MIA as they evolve from t = 0 to 1000: 
a) vorticity and b) magnetic field. 
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